%% Step 1 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Date: October 26, 2017
% Author: Produced for JdG,AL,IM by Christopher Evans at UPF 
%
% Program: imports the WIOD data for 1995 based on 2012 vintage of WIOT. 
% The data are in CSV and can be found at CSV/CSV/wiot95_row_apr12.csv.
% To avoid issues with matrix inversions, we nominally adjust 0 entries.
%
% NEW: For all code, we will change the subscripts to match those in
% master_notes3.pdf and the draft. I.e., we will change to {mn,ij} pairs,
% which represent country pair {mn} and sector pair {ij}. Global variables 
% changed and slight modification to file structure for calling up data
%
% Last Updated: 24/01/2019 by JDG
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clearvars

global alphaT rhoT lambdaT etaT rho sigma sigmaT lambda eta varphi ...
    tol maxit Pi_hat_n D_hat_n Xi_hat_mnjf Xi_hat_mnj a_hat_f a_hat ...
    N J FI_J france countrysecD firmsecD_sorted start start_sorted ...
    sum_by_country_dummy ISO ... 
    sL_n sPi_n sD_n alpha_WIOD_j alpha_WIOD_francej labour_share_PWT

%% Read data from the csv filed

% DATA: Choose which version of the wiot to import:

DATA = dlmread('CSV/WIOT_data/wiot05_row_apr12.csv',','); % WIOD 2005

% Extract data on intermediate input flows

Z = DATA(1:1435,1:1435);	

% Extract data on both flows of intemediate and final goods

WIOD = DATA(1:1435,1:1640);	

% Taking the Value Added

VA_raw = DATA(1437:1442, 1:1435); 

% Summing over the sectors so that we have VA+total intermediate
% consumption is equal to output at basic prices

VA_raw = sum(VA_raw,1);

% Final demand, gross fixced capital formation + inventories adjustment
 
PC = WIOD(1:1435, 1436:1640); 
  
% France = 15 , we aggregate BEL and LUX

france = 15;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Baseline dimensions if we also consider firms: mn,ij,f:
%
%       m: country, if placed first means exporter, 2nd is importer 
%          (production side is typically m)
%
%       n: country, if placed first means exporter, 2nd is importer 
%          (consumption side is typically n)
%
%       i: sector that exports 
%
%       j: sector that imports
%
%       f: firm. May be an exporter or importer, depending on the parameter
%          of variable we are considering (e.g., import input usage vs. 
%          export share.)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Aggregating countries (Belgium and Luxemburg), removing sectors J, P, L

% AggN: 40x41 Matrix, aggregating Luxemburg data into Belgium 

AggN= xlsread('CSV/Aggregate_Countries_Sectors.xlsx','AggC');

% AggJ: 33x35 Removing Sectors J and P from the WIOD data

AggJ= xlsread('CSV/Aggregate_Countries_Sectors.xlsx','AggS');

% N: Actual number of countries used in computation 

N=size(AggN,1);

% J: Actual number of sectors used in computation 

J=size(AggJ,1);

% NN: matrix for aggregating matrix of intermediate goods flows Z

NN=kron(AggN,AggJ); 

% FF:  matrix for aggregating matrix of final demand F

FF=kron(AggN,eye(5));  

% Z: aggregated Z matrix

Z=(NN) * (Z) * (NN');  

% F: aggregated F matrix

PC=(NN) * (PC) * (FF');   


% Set all inventories to 0:

inventory_set_zero=[1 1 1 1 0];

% PC_inventory_zero: Replicates the set to match PC, so it should be 1 for
% all final goods except inventories which is set to 0. 

PC_inventory_zero=repmat(inventory_set_zero,[N*J N]);

PC=PC.*PC_inventory_zero;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This is to calculate the alpha from the WIOD, it is done before we
%   force non-tradable sectors to have no exports and push this trade to
%   the home country 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% WIOD: Combining the intermediate and final good exports to make the WIOD
%       matrix again 

WIOD=[Z,PC]; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% R: Sums up total output from each intermediate, should equal the 
% Total Output column in excel. Needed to sum this again as we aggregated
% the data

R=sum(WIOD,2);

% VA: Value added. Recalculated because we have aggregated countries and
%     dropped sectors, therefore this will not match up exactly to the VA
%     given in the raw WIOD data.

VA=R'-sum(Z,1);

% Total_Output: Total output is the sum of all intermediate output and
%               final good output.

Total_Output= sum(WIOD,2);

% alpha_WIOD_nj: Alpha by country and by sector, this is calculated from
%                the WIOD as Value Added/Total Output

alpha_WIOD_nj=VA./Total_Output';

% Setting the NaN to zero. NaN can arise due to diving by 0 

alpha_WIOD_nj(isnan(alpha_WIOD_nj))=0;

%alpha_WIOD_nj(alpha_WIOD_nj<0)=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% WORKING WITH NON-TRADABLE SECTORS 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% List of non-tradable sectors as given in the calibration.pdf

NT_sectors={'E','F','50','51','52','H','60','61','62','63','64',...
                        'J','70','71t74','L','M','N','O','P'}';

%List of sectors in the WIOD, this is in the order of the WIOD and this is 
%the order we follow    

WIOD_sectors = {'AtB' , 'C' , '15t16' ,'17t18' , '19' , '20' , '21t22' ,...
    '23' , '24' , '25' ,'26', '27t28', '29' , '30t33' , '34t35' ,...
    '36t37', 'E' , 'F' , '50', '51', '52', 'H', '60', '61', '62', '63',...
    '64', '70', '71t74', 'M', 'N', 'O'}';                     


% Sector dummy and the sum_by_country_dummy is created to sum over country
% i and holding the sector fixed. It is a JxJ (sector) identity matrix, 
% this is repmat out to a [JxJ*N]. We use it here to sum over all sectors,
% which will be useful for our treatment of Non-Tradable sector. 

sector_dummy= eye(J);
sum_by_country_dummy= repmat(sector_dummy, [1 N]);


% Z_summed_sec: Using the dummy to hold a sector fixed and sum up over
%               export countries. Z is a [NJ NJ] and sum_by_country_dummy a
%               [J NJ} hence, the result is [J NJ]          

Z_summed_sec=sum_by_country_dummy*Z;


% Here we are looping over each of the non tradables and setting equal to 0
% the NT sector trade to foreign countries and assigning the summed by 
% sector over countries NT trade.

for m=1:N
   for i=1:J
      for j=1:size(NT_sectors,1)
         if isequal(WIOD_sectors{i},NT_sectors{j}) 
            Z(J*(m-1)+i,:)=0;                      
            Z(J*(m-1)+i,J*(m-1)+1:J*(m))=Z_summed_sec(i,J*(m-1)+1:J*(m));
         else
         end
      end
   end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% X: Combining intermediate goods and final good consumption
 
WIOD=[Z,PC]; 

% sum_of_intermediates: Sums the intermediates by n and sector j such that
%                       we have the sum of intermediates that we can use to
%                       create out gammas

sum_of_intermediates = sum(Z,1); 

% A: Input coefficients (gammas), is the intermediate output divided by the
%   sum of intermediates from country n sector j 

A=Z/diag(sum_of_intermediates+.0000000001.*(sum_of_intermediates==0));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Gamma: These are the input coefficients, taken from A (above)  

gamma=A; 

%start: Index for when we have a different country, this is because we have
%       32 sectors 
clear vars start;
start(:,1) = 1:J:size(Z,1);

% Need to add the index for the end of the matrix

start(end+1,1)= size(Z,1)+1 ;

% Working on the LHS of the equation, final goods part. PC is defined as 
% price x consumption, this is given as final goods in WIOD

PC_sectors(:,1) = 1:5:5*N;
PC_sectors(end+1)=size(PC,2)+1;

% This is summing over the 5 columns of the final good for 
% each country, and will give us PC_mnj

for s=1:N
    PC_mnj(:,s)= sum(PC(:,PC_sectors(s):PC_sectors(s+1)-1),2);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NON TRADABLE SECTORS for the final goods
%
% This sums all elements in the non tradable row and assigns it to that
% country only. It also sets the non tradable sector trade to other
% countries to zero. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Summing the final good exports by country m exported to n, used later.

PC_mj=sum(PC_mnj,2);

% Looping over each country for each sector i. This loop functions similar
% to the one before, as in, it loops over the WIOD sectors until it finds a
% non-tradable sector and then sets the whole row (all of the final good
% exports) to 0 and assigns the sum of this row (calculated previously) to
% its own country, therefore enforcing that it is non-tradable.

for m=1:N
   for i=1:J
      for j=1:size(NT_sectors,1)
         if isequal(WIOD_sectors{i},NT_sectors{j}) 
            PC_mnj(J*(m-1)+i,:)=0;
            PC_mnj(J*(m-1)+i,m)=PC_mj(J*(m-1)+i,1);
         else
         end
      end
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% PC_mj: Final goods by country m and sector j

PC_mj = squeeze(sum(PC_mnj,2));

% PC_nj: Final goods by country n and sector j 

PC_nj=sum_by_country_dummy*PC_mnj;

% PC_n: Final good consumption by country n

PC_n = sum(PC_nj,1);

% pi_c_mnj: Share of final good exports from country m to coutnry n sector
%           j. It it equal to PC_mnj/PC_nj 

PC_nj=PC_nj+0.0000000001.*(PC_nj<=0.0000000001);

for m=1:N
    pi_c_mnj(start(m):start(m+1)-1,:)= PC_mnj(start(m):start(m+1)-1,:)./PC_nj;
end

% pi_c_nj: Share of final good imports of country n sector
%           j. It it equal to PC_nj/PC_n

for j=1:J
    pi_c_nj(j,:)= PC_nj(j,:)./PC_n;
end

% Z_mnj : Intermediate trade from country m to n sector j

for n=1:N
    Z_mnj(:,n)= sum(Z(:,start(n):start(n+1)-1),2);
end

% X_mnj, total exports is intermediate output and final output 

X_mnj=Z_mnj+PC_mnj;

% X_njp: Exports from country n sector j'. This is calculated by summing
% over the importing countries k

X_nj = sum(X_mnj,2);


% X_mn: Total exports from country m to country n (any sector). This sums
%       the sectors so we have country to country.

for m=1:N
   X_mn(m,:)=sum(X_mnj(start(m):start(m+1)-1,:),1);
end

% pi_mnj: Export shares. Initially set all to 1 as this is the share of 
%          firm f in total exports from country m to country n in sector j
%          However this is the sector level equivalent to the firm
%          variable. Therefore we set export shares equal to 1 for all
%          tradable sectors. We then set export share equal to 0 for all
%          non-tradable sectors. The important exception here is that
%          non-tradable to non-tradable sector trade within the same
%          country is set equal to 1. 

pi_mnj = ones(size(X_mnj,1),N);

for m=1:N
   for i=1:J
      for j=1:size(NT_sectors,1)
         if isequal(WIOD_sectors{i},NT_sectors{j}) 
            pi_mnj(J*(m-1)+i,:)=0;
            pi_mnj(J*(m-1)+i,m)=1;
         else
         end
      end
   end
end

% Importing the one minus alpha (the labour share) from the PWT. The
% brackets [0 2 39 2] is used for the offest to get the correct columns and
% rows into Matlab. This is alpha by country, need to scale up to full
% matrix. 

labour_share_PWT = dlmread('CSV/PWTlaborshare.csv',',',[0 2 39 2]);

% countrysecD: Initialising a dummy representing each country 

countrysecD = zeros(N*J,N);

% countrysecD: Looping to put the value 1 when we country and sector equal
% to country 

for m=1:N
    countrysecD(start(m):start(m+1)-1,m)=1;
end

% alpha_WIOD_j: The alpha (labor share) taken from the WIOD (VA/Total
%               Output) and then constructed by sector using a weighted
%               average of total output

alpha_WIOD_j = sum_by_country_dummy*(alpha_WIOD_nj'.*R)./...
                        (sum_by_country_dummy*R);

% alpha_WIOD_francej: The alpha (labor share) taken from the WIOD (VA/Total
%               Output) for France

alpha_WIOD_francej = alpha_WIOD_nj(1,start(france):start(france+1)-1)';